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1.0  INTRODUCTION 


1.1  Basic  Problems 

The  plane,  steady.  Incompressible,  Navler-Stokes  equations  have  been 
subjected  to  a variety  of  numerical  attacks.  At  the  present  state  of  the 
art  there  Is  no  particular  method  that  Is  clearly  recognized  as  the  "best" 
way  to  solve  these  equations.  The  difficulties  are  varied  and  perhaps  the 
fact  that  we  cannot  easily  list  and  order  them  In  Importance  points  up  the 
largest  difficulty  of  all.  To  make  some  reasonable  start  on  studies  of 
numerical  methods  for  the  plane,  steady,  Navler-Stokes  equations  we  must 
first  classify,  at  least  In  some  rough  sense,  the  types  of  problems  that  are 
of  Interest. 

The  first  basic  classification  concerns  the  domain  (I.e.  geometry)  over 
which  the  flow  Is  of  Interest  and  the  boundary  conditions  are  to  be  Imposed.  In 
all  "external"  flows,  or  domains  extending  to  Infinity,  we  must  be  able  to 
give  proper  conditions  on  some  artificial  boundary  In  the  fluid.  In  order  to 
reduce  the  problem  to  one  In  a finite  domain.  We  do  not  attempt  to  study 
this  problem  In  our  present  work.  Instead  we  use  test  problems  In  which 
boundary  conditions  over  simple  finite  domains  are  "reasonably"  well  defined. 
For  example,  we  consider  the  driven  cavity  and  the  Inlet  region  of  a 
channel.  Even  here,  of  course,  there  are  open  questions.  But  they  do  not 
dominate  the  problem  as,  say,  the  flow  at  large  radius  does  In  the  flow  about 
a cylinder. 

The  next  basic  difficulty  concerns  the  Reynolds  number,  R»  or  rather  the 
range  of  R over  which  solutions  are  sought.  For  sufficiently  large  R 
there  will  be  thin  boundary  layers  In  the  flow,  say  of  width  6.  Thus  the 
net  spacing  and  numerical  method  must  be  able  to  resolve  variations  over 


r 
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lengths  of  the  order  of  j/10  say.  We  shall  not  attempt  to  resolve  this 
open  question  either.  But  our  methods  must  be  capable  of  working  at  large 
R If  we  can  use  appropriate  fine  nets.  Also,  we  point  out  that  most  cal- 
culations are  desired  not  for  a fixed  R value  but  rather  over  some  Inter- 
val, say  Rq  < R < Rp«  This  fact  can  play  a crucial  role  In  developing 
efficient  numerical  methods.  | 

a 

Suppose  then  that  we  have  a simple  domain  (a  rectangle  for  example)  on  | 

f 

all  of  whose  boundaries  the  "correct"  boundary  conditions  are  known.  Further,  I 

f 

steady  solutions  are  desired  for  a "moderate"  range  of  R values.  All 
numerical  attempts  to  solve  such  problems  are  Iterative.  The  literature  on 
this  subject  does  not  always  stress  these  facts  but  It  Is  clear  that  the 
dominant  concerns  are  now:  («)  to  Insure  convergence  of  whatever  Iteration 
scheme  Is  used,  (b)  to  accelerate  convergence  assuming  It  can  be  attained, 

(c)  to  Insure  that  the  converged  Iterates  accurately  approximate  a physical 
solution,  (d)  to  Insure  that  this  solution  Is  the  desired  one  (for  steady 
solutlonsmay  not  be  unique).  Of  course,  all  these  concerns  are  Interrelated 
and  may  not  be  resolved  Independently  of  each  other. 


1,^2  Basic  Methods 

\ 

In  the  present  report  and  In  our  most  recent  work  we  have  concentrated 
on  three  basic  formulations:  primitive  variables  In  -Sertlnn  ?U>fvort1c1ty- 
stream  function  UuSectlojLl^*  and  stream  functlon-biharmonic  formulation. 
In  Section  With  each  such  formulation  there  are  a variety  of  difference 
equations  that  could  be  used  and  then  there  are  a large  number  of  Iteration 
schemes  that  could  be  employed  to  solve  these  difference  equations.  The 
ultimate  goal,  of  course.  Is  to  find  the  best  combination  of  all  these 
techniques. 


2 


We  have  discontinued,  for  the  present,  our  earlier  primitive  variable- 
splitting  technique  studies^^^.  They  have  some  very  attractive  features  but, 
as  we  had  used  them,  they  seemed  quite  sensitive  to  the  boundary  treatment. 
Also,  their  convergence  rates  did  not  compare  to  that  of  our  new  b1 harmonic 
method . 

The  primitive  variable  methods  of  SectionZ.O  are  based  on  the  reportedly 
successful  techniques  of  Spalding^^^  and  his  coworkers.  Basically,  ADI 
methods  are  used  to  compute  the  velocity  from  the  momentum  equations  and  then 
these  velocities  and  the  pressure  are  adjusted  to  satisfy  continuity.  This 
latter  procedure  is  closely  equivalent  to  solving  a Poisson  equation  for  the 
pressure.  This  is  also  done  by  an  inner  loop  of  ADI.  Careful  tuning  of 
this  method  seems  to  be  required  and  adjustments  of  various  kinds  are  made 
when  it  is  applied  to  different  problems.  It  is  not  clear  that  any  analysis 
of  convergence  rates  or  order  of  accuracy  can  be  done  for  this  complicated 
method. 

The  vorti city-stream  function  formulation  of  Section  3*0  is  classical 
for  plane  steady  viscous  flows.  With  more  recent  developments  in  numerical 
analysis.  It  should  be  possible  to  make  this  formulation  much  more  efficient 
than  it  has  been  In  the  past.  Thus  we  use  ADI  on  the  vorticity  equation, 
suitably  linearized,  and  SOR  on  the  Poisson  equation  (vorticity  definition). 
There  are  the  usual  troubles  with  boundary  conditions  in  this  method  since 
vorticity  is  not  known  on  the  boundaries  and  is  In  fact  generated  there.  It 
seems  clear  that  the  use  of  fast  Poisson  solvers  could  greatly  enhance  the 
efficiency  of  these  procedures;  we  hope  to  Investigate  this  In  the  future. 

Finally,  In  Sect1on4.0  we  use  the  biharmonic  or  fourth-order  formula- 
tion In  terms  of  a stream-function  (i.e.  eliminate  the  vorticity).  Now 
there  are  no  problems  with  the  proper  boundary  conditions  (on  rigid  walls. 


for  Instance).  Also*  there  are  not  coupled  systems  of  differential  equations 
to  be  solved.  Newton's  method  and  modifications  of  It*  In  particular  the 
Newton-Chord  method*  are  extremely  fast  for  solving  the  finite  difference 
equivalents  of  this  formulation.  The  main  drawbacks  with  this  method  are 
the  possibly  large  storage  requirements  and  the  time  It  takes  to  get  the 
LU -decomposition  of  the  large  band-structured  coefficient  matrix  that  occurs. 
But  tentative  tests  Indicate  that  this  scheme  can  be  as  fast  as  the 
others  and  works  for  a far  greater  range  of  Reynolds  numbers. 
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2.0  PRIMITIVE  VARIABLE  FORMULATION 
2.1  Equations  and  Grid 

Of  the  few  methods  which  solve  the  Incompressible  Navler-Stokes  equations 
In  primitive  variables,  the  most  successful  one  derives  from  the  work  of 
Chorln^^^,  Amsden  and  Harlow^^^,  and  Spaldlng^^^.  This  entire  group  can  be 
described  as  being  “pressure  corrector"  methods.  The  meaning  of  this  phrase 
will  become  clear  In  the  detailed  description  which  follows. 

The  most  familiar  of  these  methods  are  the  ones  which  have  been  developed 
at  Imperial  College.  The  latest  in  the  series  of  computer  programs  which 
utilize  this  pressure  corrector  scheme  Is  known  as  TEACH.  Now,  although  the 
basic  method  is  the  pressure  corrector  approach,  numerous  computational  devia- 
tions from  a straightforward  application  of  the  method  have  been  incorporated 
into  the  TEACH  code  to  render  it  useful  for  solving  practical  problems.  These 
devices  range  from  the  incorporation  of  upwind  differencing,  and  different 
under-relaxation  parameters  for  different  variables,  to  the  Inclusion  of 
"false  source"  terms  In  the  equatlons^®^  In  order  to  enhance  convergence,  and 
the  use  of  a block  correction  procedure^®^  to  correct  the  difference  equations 
for  discretization  errors.  All  these  tend  to  obscure  the  basic  properties  of 
the  pressure  corrector  method  itself  and  hence  make  comparisons  with  other 
methods  difficult,  so  the  procedure  described  below  is  the  basic  concept  of  the 
pressure  corrector  method,  applied  to  various  test  problems,  but  bereft  of  all 
computational  devices. 

Previous  Investigations  have  provided  related  Information  concerning 
Implementation  of  the  method  (see  Spaldlng^^^  or  Roache^^^  for  example)  but  It 
Is  useful  to  have  a complete,  detailed,  development  In  one  self-contained  work, 
and  much  of  the  interpretation  and  some  details  which  are  Included  here  are  new. 

The  equations  to  be  solved  are  the  suitability  nondimenslonal  Incompressible 
form  of  the  Navler-Stokes  equations,  I.e. 
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3U  3U  . dU 

•~  + U — + V — 

3t  3X  3y 


IV  + y 3V  y ^ 
3t  3X  3y 


1 / 3^U  . 3^U  \ 

n;?  V) 

1 /3^V  ^ 3^v\ 


(2.1) 


(2.2) 


1“+  3V  = 0 

3X  3y 


(2.3) 


In  these  equations  the  time  derivative  Is  regarded  as  a means  of  approaching 
the  asymptotic  steady  state,  and  a consistent  time  development  of  the  equations 
will  not  be  considered.  The  TEACH  method  considers  the  steady  equations,  and 
solves  this  elliptic  system  by  relaxation  methods.  The  conservation  form  of 
these  equations  can  be  shown  to  possess  the  conservative  property^^^,  hence 
conserve  mass  In  the  global  sense  when  cast  In  finite-difference  form,  so 
Instead  of  eqs.  (2.1)  and  (2.2),  the  following  equations  are  used: 


If  * k I7 


(2.4) 


(2.5) 


The  finite-difference  form  of  these  equations  Is  solved  on  the  staggered 
grid  depicted  In  fig.  1.  Arguments  as  toihe  physical  basis  for  this  grid  have 
been  made,  and  It  has  become  standard  practice  for  all  primitive  variable 
solution  techniques  (Including  finite-element  methods^^^).  However,  It  should 
be  noted  that  the  original  pressure  corrector  method  of  Chorin  did  not  use 
this  grid,  but  the  standard  one  with  all  variables  defined  at  the  same  point. 

Note  the  location  of  the  boundaries  of  the  solution  domain  and  the  posi- 
tion of  each  variable  relative  to  the  grid  and  boundaries.  On  this  grid 
the  Interior  values  of  the  pressure,  p^  j,  are  defined  on  2 ^ 1 ^ IMAX 
and  2 £ j £ JMAX.  The  values  along  grid  lines  1 or  MAX+1  must  be 
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determined  by  boundary  conditions,  or.^as  will  be  seen  later,  other  means. 

The  Interior  values  of  the  u-velocity  are  defined  on  3 <_  1 ^ IMAX,  2 J JRAX. 

The  values  at  1=2  and  1 = IMAX+1  fall  directly  on  the  boundary  lines  of 
the  solution  domain  and  are  provided  by  boundary  conditions.  The  values  at 
j = 1 and  j - JMAX+1  are  exterior  to  the  solution  domain  and  are  specified 
In  combination  with  Interior  values  of  u to  provide  the  correct  boundary  con- 
dition at  j = 3/2  and  j » JMAX  + 1/2.  The  values  at  1 = 1 are  never  used. 
The  situation  for  the  v-velocity  Is  Identical  to  u with  the  roles  of  1 and 
j Interchanged.  Thus  v Interior  values  are  2^1  ^ IMAX,  3 j ^JMAX; 
j “ 2 and  j = dMAX+1  values  fall  on  boundaries.  1 » 1 and  1 ■ IMAX+1  are 
exterior  and  calculated  In  combination  with  Interior  values,  and  j ■ 1 values 
are  never  used.  The  exact  specification  of  the  boundary  values  will  be  given 
below  In  section  2,4, 

2.2  Finite-Difference  Equations 

The  spatial  differencing  scheme  to  be  used  will  be  central  differencing 
at  all  points;  this  Is  standard  practice.  However,  we  choose  an  Implicit 
Integration  technique,  as  In  the  formulation  for  vorticity  streamfunctlon  vari- 
ables, on  the  assumption  that  rather  arbitrary  choices  of  time  steps  may  be 
taken  to  enhance  the  convergence  to  the  desired  steady-state  solution.  This 
point  will  be  returned  to  later.  The  ultimate  Integration  procedure  will  be 
an  alternating  direction  Implicit  (ADI)  method,  but  due  to  the  nature  of  the 
way  the  pressure  corrector  equations  are  derived,  the  ADI  method  will  be 
arrived  at  In  an  Indirect  way.  The  main  reason  this  Is  necessary  Is  that  the 
pressure  corrector  equations  are  derived  from  considering  the  1 1nearl  zed 
difference  equations  and  not  the  original  differential  equations. 
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To  derive  the  final  form  of  the  equations,  consider  the  fully-implicit 
differencing  of  equations  (2.4)  and  (2.5).  The  x-momentum  equation  is 
centered  at  the  point  (i  -1/2,  j),  and  the  y-momentum  equation  is 
centered  at  (i  • j -1/2).  The  two  equations  are 


Ay 


= _!lJ *^<-1.3  . 1 n*1/2.J  J ^ “)-3/2..1 

4X  R \ 


AX 

+ 1-1/2J^T  ^“i-l/2.j;^  *^i-l/2.i-l 


) 


(2.6) 

2\n+l  / ,2,n+l 


- . c">r^R.3.i/R  - , (v^)ri  - 

(At/Z)  ax  ay 


ax 

n+1  _ n+1  / n+1  _ . n+1  n+1 

■ -"u  ’’l.J-l  . 1 |Vl. 3-1/2  ^*1 .1-1/2  * ''l-l .1-1/2 

iy  r\ 


^ .3*1/2  ^*1 .3-1^2  * ''l.J-3/2 


4y 


) 


(2.7) 


where  (n+1)  is  the  new  time  station,  (n)  is  the  old  time  station,  the  factor 
of  1/2  in  the  time  step  is  for  convenience  later,  and  the  pressure  for  the  moment 
is  assumed  to  be  known.  After  multiplying  through  by  at  we  can  define  the 
following  quantities. 


r “ At  P _ At 

''X  " AX  * V " Ay  ’ 


D * , 

^ Rax'^ 


At 

5 

RAy^ 


(2.8) 


Since  the  convective  terms  in  eqs.  (2.6)  and  (2.7)  are  nonlinear  and  we  can 
only  solve  linear  difference  equations,  we  must  linearize  the  equations. 

The  linearized  values  of  the  dependent  variables  will  be  denoted  by  an 
overbar.  Dealing  with  only  eq.  (2.6)  for  the  present,  we  get 
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ii 


i! 


* S^''i-l/2.j+l/2‘'i-l/2,j+l/2  ''i-l/2.j-l/2“i-l/2.j-l/2^ 


+ Dy(ujtl/2,j+l 


^“1-1/2,^  ^ “i-l/2.j-l^ 


(2.9) 


Since  the  values  of  u are  only  known  at  the  half-integer  values  of  i,  and 
the  values  of  v are  known  only  at  the  integer  values  of  j on  the  staggered 
grid,  simple  averages  (which  are  formally  second-order  accurate)  are  used  to 
eliminate  these  values  in  eq.  (2.9), 


2 fi-l/2,j  ”“i-l/2,j)'^  V“i,j  ' 2 ^“i-l/2,j  ^ “i*l/2,j^ 


“i-Ko  * \ ^“i-l/2,j  “i-3/2,j^^ 


. - r-  1 /„n+l  . „n+l  N 

V i-l/2J+l/2  * I '“i-l/2,j+l  “i-l/2,j^ 


''i-l/2,j-l/2  ' I ^“i-l/2,j  * ‘'i-l/2,j-l^^  “ ‘^x^Pi,j 


r,n+l  \ 

Pi-1  ,y 


Px'“i+l/2,j  ^‘'i-l/2,j  i-3/2,j^ 


'^y^'^i-l/2,j+l  ^“i-l/2,j  ‘'i-l/2,j-l^ 


(2.10) 


where  the  averages  are  used  in  calculating  the  linearized  (barred)  quantities 
also. 

After  regrouping  the  terms  in  eq.  (2.10),  the  final  form  of  this  differ- 
ence equation  can  be  written 
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^N“i-1/2J+1  * ^‘'i-3/2,j  ^2  + Apj^  + ApY)uJ_y2,j  * ^E“l+l/2,j 


where 


* ?S''l-l/2J+V2  °y»  “ r Vl-l/2,j-l/2  “‘'y 


2Vi-l.j  °x’ 


(2.12) 


iU  ,1 


^PX  * 2 ^x^“i,j  “i-1.j^  ^'^x’  ^PY  “ 2 ^y^''i-l/2,j+l/2  ''l-l/2J-l/2^ 


Similarly,  the  final  form  of  the  v*momentum  equation,  (2.7),  can  be  written 


.V..n+1 


iV.,n+l 


A*v""'  + + f2  + A*  + a'^ 

Vl,j+l/2  Vi-l,j.l/2  ' ''PX  ''PY^  1,j-l/2  Vl+l,j-l/2 


iV..n+l 


^ Vi  ,j-3/2  ''l  ,j-l/2  ^ ,j  Pi  ,j-l  > ° 


n+1  _n+l 


where 


(2.13) 


aJJ  « -y  Vi -1/2, j -1/2  “°x’ 


.V  _ 1 


^S*  rVu-l  "^y 
“ ? Vl+l/2,j-l/2  "‘^x 


(2.14) 


^PX  “ 2 ^x^“l+l/2,j-l/2  “l-l/2,j-l/2^  * 


ApY  - 5- Cy(Vi  ,j  ''^,j.l)  + 2Dy 


The  only  remaining  equation  the  continuity  equation,  (2.3).  This  Is 
placed  In  finite-difference  form  by  centering  the  equation  at  the  point  (1,j). 


n+1  _ n+1  + ^ / n+1  _ n+1  x ^ q 

“l+l/2,j  “l-l/2,j  * Ay  '''l,j+l/2  ''l,j-l/2^  ^ 


(2.15) 
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These  three  equations,  (2.11),  (2.13)  and  (2.13)  form  the  basis  for  the 
pressure  corrector  method. 

The  Implementation  of  the  method  follows  from  the  following  argument. 

If  the  correct  pressure  were  known,  then  eqs.  (2.11)  and  (2.13)  could  be 
solved  for  u and  v,  and  eq.  (2.15)  would  be  satisfied  exactly.  However, 

It  Is  not  known,  so  we  have  to  make  a guess  for  It,  and  then  correct  It  In 
a rational  way  to  satisfy  eq.  (2.15)  which  we  want  to  hold.  Hence,  we  assume 
that  the  true  solution  Is  given  by 


u^^  » u*  + u' 


= V*  + v' 


(2.16) 


_n+1  _ A . I 

p = p»  + p 

where  the  starred  quantities  are  the  guesses  and  the  primed  quantities  are  the 
corrections.  Placing  (2.16)  Into  eq.  (2.11)  for  u yields 

* « * ''PX  * *PY><“*  * "')l-l/2.j 
* *e'“*  * “'•ltl/2,j  * ''s*"*  * “''l-1/2J-l  ““?-l/2,j 

+ C,[(p*  ♦ - (p*  + P')i.,,j]  = 0 

Hence,  the  equation  for  the  predicted  (guessed)  value  of  u*  Is  Identical  to 
equation  (2.11)  but  with  (n>1)  values  replaced  by  starred  values,  since 
eq.  (2.17)  Is  linear  and  can  be  split  Into  a predictor  equation  and  a corrector 
equation.  The  remaining  terms  In  eq.  (2.17)  give  the  full  equation  for  the 
correction  value  u': 
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PY 


)u' 
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''E‘'i+l/2.j  * ^S“i-l/2.j-l  J “Pji  • 0 (2.18) 

This  equation,  as  It  stands.  Is  of  no  value,  but  an  approximate  form  of  the 
equation  leads  to  the  final  solution.  For  more  details  of  the  following 
argument,  see  Chorln^^^  or  Amsden  and  Harlow^^^.  In  order  to  alter  the 
velocity,  but  not  change  the  value  of  the  vorticity,  the  velocity  can  only 
be  corrected  by  the  gradient  of  a scalar  function.  Hence,  eq.  (2.18)  Is 
approximated  by 


(2.19) 


In  an  analogous  way  for  the  v-momentum  equation.  It  can  be  shown  that  the 
predicted  value  of  v*  Is  obtained  from  eq.  (2.13)  with  (n-i-l)  replaced  by 
starred  values,  and  the  correction  for  v front 


''l'.J-l/2 


(2  + A 


PX  ^PY' 


(2.20) 


Note  In  passing  here  that  had  the  values  for  a predicted  and  corrected  value 
been  placed  In  either  the  full  differential  equations  (2.4)  and  (2.5)  or  the 
nonlinear  difference  equations  (2.6)  and  (2.7)  the  results  for  the  correction 
values  would  have  been  significantly  more  complicated. 

To  determine  the  equation  for  the  correction  (primed  values), 
account  for  the  fact  that  the  continuity  equation  has  not  yet  been  satisfied. 
Thus,  placing  eqs.  (2.16)  for  u”^^  and  v*^^^  Into  the  differenced  continuity 
equation,  (2.15),  and  using  the  relations  (2.19)  and  (2.20)  for  the  correction 
velocity  we  get 
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Ay  Vi  j+i  Vi-i  j * ®p  * ^ * ®p^^Pi  ,j  ~ ®EPi+i  ,j 


(2.21) 


where 


(2  + aJj(  + AJy) 


.(2  + aJjj  + aJy) 


(2  + aJx  + A^ 


(2  + aJx  + AJy) 


(2.22) 


“ “T+l/2.j  “l-l/2.j  ^ f7  ^''i  ,j+l/2  “ ''T.j-1/2^ 


Clearly  (2.21)  Is  a difference  form  of  a Poisson  equation  for  the  pressure. 


2.3  Solution  Algorithm 


The  solution  procedure  for  this  predictor-corrector  scheme  is  as  follows: 
(a)  Guess  a value  of  the  pressure,  e.g.  the  pressure  at  the  last  integration 
station;  (b)  compute  the  linearized  quantities  present;  also  from  last  inte- 
gration station;  (c)  integrate  the  two  momentum  equations,  subject  to  the 
imposed  boundary  conditions  discussed  below  in  section  2.4,  to  get  the  values  of 
u*  and  V*;  (d)  find  p'  from  the  Poisson  equation,  (2.21);  (e)  correct 

the  starred  quantities  using  eq.  (2.16)  to  get  the  new  values  of  u,  v,  p; 

(f)  continue  this  process  until  a satisfactory  convergence  is  obtained. 

Step  (c)  above  requires  the  solution  of  the  two  momentum  equations  which 


now  have  the  form 


+ aHu? 


+ (2+  A^„  + A^y)u* 


Vl-l/2,j+l  ''W“i-3/2.j 


+ aSu*. 


PX  "PY'“i-l/2,j  "E“i+l/2,j 


^ ^S“i-l/2,j-l  ““l-l/2,j  * ^x^Pi,j  "Pi-l,j^  “ ° 
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Vl,j+l/2  * Vl-lJ-l/2  * ^ 


PX 


V^''i,J-l/2  * Vl+lJ-l/2 


* ^s''lJ-3/2  ''i,j-l/2 


(2.24) 


Rather  than  solving  these  in  this  fully  implicit  form.in  which  they  were 
written  for  convenience  in  deriving  the  predictor  and  corrector  steps  of 
the  methods*  an  alternating  direction  procedure  will  be  used  to  generate  a 
scheme  where  only  tri diagonal  matrices  need  to  be  solved  for  the  unknown 
values.  Writing  eqs.  C2.23).  and  (2.24) in  a more  convenient  shorthand 
notation 

+ 0.0  (2.25) 


where  the  correspondence  between  (2.25)  and  either  (2.23)  or  (2.24)  can  be 
found  by  inspection,  and  only  the  relative  (1,j)  positions  of  the  unknowns 
on  the  difference  mesh  of  figure  1 has  been  shown.  The  factor  2,  in  the 
middle  term  above,  which  was  Introduced  into  the  time  differencing  of  eqns. 
(2.6)  and  (2.7)  will  now  be  seen  to  allow  a symmetric  AOI  scheme  to  be  formed. 
Eq.  (2.25)  can  be  split  into  the  following  ADI  scheme 

(1  + Fx)q  + (-1  + Fy)q"  + D = 0 (2.26a) 

(1  + Fy)  ^ (-1  + Fj^)q  + D - 0 (2.26b) 

which  is  equivalent,  to  second-order  accuracy,  to  a Crank-Nicolson  scheme 
expressed  as 

q -q"  + (Fx  + Fy)(q  + q")  + 2D  - 0 


(2.26c) 


where 

fx  ' ‘•x’l-l  * Vi  * Vi+1 

' Vm  * Vj  * Vj+1 

Thus,  equations  (2.26)  can  be  looked  upon  as  a Peaceman-Rachford^*'^  AOI 
solution  procedure  for  either  u*  or  v*,  and  this  only  requires  the  use  of 
a two-sweep,  or  Thomas,  algorithm  for  the  Inversion  of  tridlagonal  matrices. 

The  Poisson  equation  for  the  pressure  correction,  eq.  (2.21),  needs  to 
be  solved  In  step  (d)  above.  In  the  time  asymptotic  solution  of  steady-state 
problems.  It  has  been  found  that  the  Poisson  equation  need  not  be  solved 
exactly  at  each  fictitious  time  step,  but  that  only  a reasonably  consistent 
value  of  p'  need  be  generated  by  the  Poisson  solver.  This  Is  accomplished 
by  turning  the  elliptic  Poisson  equation  Into  a nonhomogeneous , parabolic 
equation  for  p',  and  advancing  this  equation  a few  (usually  3-5)  time  steps 
for  the  value  of  p'.  Hence  the  Poisson  equation,  (2.21),  can  be  Immediately 
cast  Into  a form  exactly  equivalent  to  eq.  (2.26)  by  adding  2/&x  p^  ^ and 
subtracting  2/ At  p^  ^ to  equation  (2.21).  Now  the  ADI  formalism  can  be 
used  to  march  the  converted  eq.  (2.21)  three  to  five  steps  for  a consistent 
p'  and  the  solution  procedure  outlined  at  the  beginning  of  this  section  can 
proceed. 

2.4  Boundary  Conditions 

Although  the  use  of  the  primitive  variable  formulation  allows  the  direct 
specification  of  the  velocity  values  on  the  boundaries,  rather  than  indirectly 
as  In  the  case  of  the  vorticity  and  stream- function  formulation,  the  use  of 
the  staggered  grid  depicted  In  figure  1 complicates  things.  Were  all  velocity 
points  coincident  with  the  boundary, specification  would  be  trivial.  However 


(2.27a) 

(2.27b) 


I 


since  half  the  velocity  nodes  needed  to  specify  boundary  conditions  lie 
outside  the  computational  domain,  special  relations  are  necessary  to  provide 
the  correct  boundary  values  at  the  midpoint  between  the  velocity  nodes.  On 
the  other  hand,  the  boundary  conditions  on  the  pressure  (actually  the  pres- 
sure correction,  as  will  be  seen)  are  quite  simple,  and  easily  applied. 

In  this  section,  the  four  boundaries  of  the  domain  shown  In  figure  1 
are  denoted  by  the  obvious  terms  upper,  lower,  left  and  right  boundaries. 
Conditions  for  u on  each  of  these  will  be  described  first,  followed  by 
conditions  for  v on  each.  Next  will  be  a discussion  of  the  pressure 
boundary  conditions,  and  the  consequences  on  the  pressure  of  using  the  pres- 
sure corrector  method.  Finally,  the  Incorporation  of  typical  boundary  condi- 
tions Into  the  Implicit  solution  algorithm  described  In  section  2.3  will  be 
shown. 

The  u-velocity  nodes  fall  directly  on  the  boundaries  on  the  left  and 
right  sides  of  the  domain.  Thus,  If  the  velocity  Is  given  on  these  boundaries, 
for  example  a solid  wall,  or  a specified  Inflow  profile,  the  value  Is  fixed 
at  the  given  value.  If  outflow  Is  specified  at  the  right  boundary,  I.e. 
au/dx  » 0,  then  the  following  relations  are  used  to  connect  the  variation 
of  the  boundary  velocity  with  two  Internal  points. 

“lMAX+l,j  “ ^^“lMAX,J  “ “lMAX-l,j^^^ 

This  relation  Is  second-order  accurate  and  does  not  destroy  the  tridlagonal 
system  of  equations  when  Implemented  as  shown  below. 

For  all  problems  discussed  here,  the  upper  and  lower  boundaries  can  be 
considered  solid  walls.  Thus  the  u velocity  on  these  walls  Is  given  by 
the  wall  velocity  due  to  the  no-sllp  condition.  Unfortunately,  the  boundary 
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Is  located  between  u nodes,  see  figure  1,  and  hence  the  u-value  given  there 
cannot  be  directly  specified.  Thus,  a relation  giving  u at  the  boundary 
In  terms  of  the  exterior  and  Interior  points  must  be  derived.  The  standard 
relation  which  has  been  recommended  for  this  problem  Is  based  on  a simple 
average*  For  the  top  wall,  denoting  the  velocity  at  the  wall  by  subscript 
w,  we  can  write 


T ^“l ,JMAX+1 


* “l,JMAX^ 


Hence  the  boundary  condition  to  be  used  for  the  unknown  velocity  at  node 
JMAX+I  can  be  expressed  as 

“l ,JMAX+1  “ 2“w  “ “l ,JMAX 

This  Is  a formally  second-order  accurate  result.  However,  a higher-order 
relation  which  still  retains  the  tridlagonal  form  can  be  derived.  This 
gives 


“l ,JMAX+1  “ I “w  ^‘^1,JMAX  3 “l,JMAX-l 

Now,  althou'  ..)e  recommended  relation,  (2.29)  1$  formally  second-order 
accurate,  consider  the  case  of  a fully  developed  pipe  (Polseullle)  flow. 

Here  the  u velocity  has  a parabolic  form.  Equations  (2.29)  Implies  that 
the  velocity  variation  around  the  wall  point  Is  linear,  which  In  the  case 
of  Polseullle  flow.  It  certainly  Is  not.  Relation  (2.30)  was  derived  by 
allowing  a quadratic  variation  around  the  wall  and  Is  consistent  with  a para- 
bolic profile,  such  as  Polseullle  flow.  The  use  of  these  on  the  test  prob- 
lems will  be  discussed  In  section  2.5  where  the  results  are  presented. 

The  situation  for  v boundary  conditions  Is  very  similar  to  what  was 
just  stated  for  u above,  except  with  1 and  J reversed.  The  upper  and 
lower  walls  coniclde  with  v nodes  ar.ti  so  v can  be  specified  here  exactly. 
For  all  cases  considered  here  v * 0 on  both  walls.  The  left  boundary  Is 
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either  a wall  or  a specified  Inflow,  and  so  Is  known.  Thus  an  appropriate 
version  of  eqs.  (2.29)  or  (2.30)  for  v at  1 » 1 can  be  used.  The  right 
boundary.  If  It  Is  a wall,  uses  eq.  (2.29)  or  eq.  (2.30)  with  u replaced  by 
v and  1 and  J unchanged.  If  It  Is  an  outflow  boundary,  then  3v/3x  - 0 
must  be  used  In  order  to  specify  j.  Here,  the  standard  relation  given 

Is  the  only  one  which  could  be  found  to  satisfy  this  gradient  condition,  con- 
sistent with  maintaining  the  tridlagonal  form  of  the  system,  I.e., 

''iMAX+l.j  * ''iMAX.J 


q 

} 

i 


i 
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The  boundary  conditions  necessary  for  the  solution  of  the  pressure  cor- 
rector equation,  (2.21)  follow  directly  from  the  definitions  of  the  velocity 
corrections,  eqs.  (2.19)  and  (2.20).  Since  In  all  cases  here  the  velocity 
Is  to  be  specified  on  every  boundary,  we  solve  the  equations  for  the  predicted 
(starred)  velocities  by  setting  u*  and  v*  equal  to  the  correct  boundary 
conditions.  Hence,  the  correction  (primed)  velocities  are  all  zero  on  the 
boundaries.  However,  using  eqs.  (2.19)  and  (2.20)  these  correction  velocities 
are  related  to  the  gradient  of  the  correction  pressure,  so  at  the  left  and 
right  boundaries  where  u Is  specified,  we  have  dp'/sx  « 0,  and  on  the  upper 
and  lower  boundaries  dp'/sy  * 0.  Thus  the  Poisson  equation  for  the  pressure 
correction  Is  solved  subject  to  Neumann  boundary  conditions. 

Note  that  these  conditions  are  on  the  pressure  correction,  p',  and  not 
the  guessed  p*  or  corrected  pressure  p.  In  fact.  In  the  formulation  used 
and  described  In  this  report,  the  boundary  values  of  p or  p*  are  never 
used.  Checking  the  differencing  and  the  grid  will  cdnfirm  this.  This  fact, 
which  does  not  seem  to  have  been  emphasized  before.  Is  a plausible  result  of 
solving  the  equations  In  primitive  variables.  Most  pressure  boundary  condi- 
tions are  derived  from  the  equations  of  motion  since  the  actual  pressure  on 
a boundary  Is  not  known  a priori  In  any  flow  situation.  They  are  considered 
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"extra"  conditions  necessary  to  complete  the  particular  mathematical  (numer- 
ical) formulation  being  used.  In  fact,  they  are  extra,  because  In  the  vorticity 
stream- function  formulation  a fourth-order  biharmonic  equation  for  the  stream 
function  can  be  generated.  This  requires  only  two  boundary  conditions  at  each 
boundary  surface.  They  are  the  function  and  normal  derivative,  which  are  In 
fact  u and  v.  No  further  conditions  are  necessary.  Essentially  the  same 
thing  Is  done  In  the  "pressure  corrector"  method.  The  velocities  are  given 
as  conditions,  and  then  a consistent  set  of  boundary  conditions  on  a non- 
physical variable  Is  derived  and  used  to  compute  the  relevant  physical  quantity, 
I.e.,  the  pressure.  Never  Is  any  recourse  needed  for  a "pressure"  boundary 
condition.  In  this  way,  the  pressure  corrector  scheme  models  the  physics  of 
the  flow  closely  In  that  an  Internally  generated  value  of  the  pressure  is 
obtained  at  all  points  Inside  the  solutipn  domain  without  any  reference  to 
boundary  values  of  pressure.  If  the  value  of  the  pressure  on  the  boundary 
(which  lies  between  pressure  nodes)  Is  desired,  an  extrapolation  from  Interior 
points  can  be  used. 

To  Incorporate  any  of  these  boundary  conditions  Into  the  solution  algo- 
rithm the  following  method  Is  used,  which  keeps  the  tri diagonal  form  of  the 
matrix  Intact,  and  allows  the  use  of  the  tridiagonal  solver  given  In  the 
Appendix.  Using  the  notation  of  the  Appendix,  the  general  form  of  the 
difference  equations  generated  In  this  section  are 

\^-i  ■*”  * ^kVi  ° 

where  the  Q are  the  unknown  variables  and  the  A,  B,  C,  D are  coefficients 
computed  at  the  K location.  It  Is  presumed  that  boundary  values  for  Q 
(either  computed  or  prescribed)  are  given  at  a starting  location,  KST,  and 
an  ending  location,  KND.  The  tridiagonal  solver  Is  then  Invoked  from  KST+1 
to  KNO-1.  Thus,  at  the  beginning  and  final  stations  of  the  solver's  use. 
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eq.  (2.32),  appears 

^KST+l^KST  * ®KST+1^KST+1  * ^KST+l^KST+2  * “^KST+l 


* “ksw  ' “ 

(2.33a) 

* '^KND-1  “ ° 

(2.33b) 

However,  at  each  of  these  stations,  the  boundary  value  1s  either  known  out- 
right, or  Is  given  In  terms  of  the  two  next  Interior  points.  As  an  example 
writing  boundary  condition  (2.30)  In  general  form  consistent  with  (2.32)  and 
«'2.33)  gives 

\nD  “ '^‘^KNO-1  ^ *^KN0-2  * ^ 

This  Is  In  fact,  the  most  general  boundary  condition  which  can  be  Incor- 
porated, including  all  conditions  from  fixed  values,  for  which  b • c = 0,  to 
forms  like  eq.  (2.30).  Placing  eq.  (2.34)  Into  (2.33b)  gives  the  following 
equation  at  the  last  point  where  the  solver  Is  employed. 


where 


^KND-l^KND-2  * ®KND-1^KND-1  ^ \nD-1  “ ° 


^KND-1  ' ^KND-1  * ®^KND-1 

®KND-1  * ®KND-1  ^ “^^KNO-l 

°KND-1  “ °KND-1  * ‘^^KND-l 


(2.35a) 

(2.35b) 

(2.35c) 

(2.35d) 


Thgs  eq.  (2.35a)  Is  In  the  correct  form  for  the  solver  to  proceed,  and  there 
Is  a similar  equation  at  KST+1.  All  that  Is  required  Is  that  the  coefficients 
A,  B,  C,  D be  modified  to  Incorporate  the  boundary  values  as  In  (2.35b)  - 
(2.35c)  at  the  two  end  points  KST+1  and  KND-1 , and  after  the  solver  gener- 
ates the  values  of  Qj^  (KST+1)  £ K ^ (KND-1),  a separate  calculation  gener- 
ates and  from  the  definitions  of  the  boundary  values,  such  as 

eq.  (2.34). 
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3.0  VORTICITY  STREAM  FUNCTION  FORMULATION 


3.1  Equations  and  Grid 

The  equations  used  to  solve  for  the  vorticity  and  stream  function 
as  dependent  variables  are  easily  derived  from  the  Incompressible  primitive 
variable  form  of  the  Navler  Stokes  equations.  The  vorticity,  for  the  two 
dimensional  case.  Is  defined  by  the  single  component 


3U  { 

^ " ay  ax  ' 

An  equation  for  this  variable  can  be  obtained  by  cross-differentiating 
equations  (2.1)  and  (2.2)  and  eliminating  the  pressure  to  yield 


(3.1) 


“ IJ-  ^ ly  “ 1 (l^  0) 


(3.2) 


Again,  the  time  derivative  Is  Included  In  the  equations  as  a means  of 
reaching  the  asymptotic  steady  state  which  Is  desired;  no  true  time  dependent 
solutions  are  sought.  Using  the  definition  of  the  stream  function  In  two 
dimensions 


u.|t  . v.-|t  (3.3) 

the  continuity  equation,  eqn  (2.3)  Is  automatically  satisfied,  and  the 
equation  which  determines  t|»  can  be  found  by  placing  eqn  (3.3)  Into 
definition  (3.1)  to  yield 


(3.4) 
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These  tMO  equations  (3.2)  and  (3.4),  along  with  the  auxiliary  relations 
(3.3)  for  the  velocity,  are  the  equations  to  be  used  In  this  formulation 
of  the  Incompressible  Navi er- Stokes  equations. 
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The  finite-difference  form  of  these  equations  Is  solved  on  a standard, 
five  point,  mesh,  depicted  below 

(1J+1) 

(1-l.j)  tlJ)  (1+lJ) 

• • • 

(U-l) 

In  contrast  to  the  primitive  variable  procedure,  here  all  the  dependent 
variables  are  located  at  the  same  points,  and  the  boundaries  of  the  solution 
domain  pass  through  the  first  and  last  row  or  column  of  nodes.  Thus  the 
boundary  values  are  Imposed  at  the  boundary,  and  no  Interpolation  Is  required. 

However,  the  boundary  conditions  for  the  vortlclty  stream  function 
formulation  still  are  not  applicable  In  a straightforward  manner.  We  need 
conditions  on  all  boundaries  for  ^ and  c.  What  we  have  are  no  slip  condi- 
tions for  u and  v.  The  normal  velocity  condition  can  be  translated  Into 
a condition  for  »|»,  but  the  vortlclty  boundary  condition  Is  still  lacking. 

The  method  for  determining  the  boundary  vortlclty  Is  described  In  Section  3.4 
below,  and  Is  one  of  the  standard  means  of  accomplishing  this. 

3.2  Finite  Difference  Equations 

Both  the  ’time  dependent'  equation  for  vortlclty  transport,  (3.2),  and 
the  Poisson  equation  for  the  stream  function  (3.4)  are  solved  by  using 
central  differences  In  space.  Unlike  the  primitive  variable  solution 
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procedure,  these  equations  are  solved  directly  for  the  dependent  field 
variables,  and  so  the  time  Integration  In  eq.  (3.2)  can  be  done  Immediately 
by  an  ADI  procedure.  Since  the  convection  turns  on  the  left  hand  side  are 
non-linear  they  are  handled  In  much  the  same  way  as  In  the  primitive  variable 
section,  where  an  overbar  was  used  to  denote  a dependent  variable  which  has 
been  linearized.  Thus,  the  difference  equations  appear: 

r*  . — rP  . /r*  _ V . n •> 


, - /'1+i.j  "'1-i.j  \ 7 /ci;  ..,1  -cj . i\ 

At/2  + “l,j\  * ''l,j 


(3.5a) 


W j 


4P 


n+1  ,n+l  \ 


(3.5b) 


- 2«*,j  ^ «1-1 


MJ+l 


.n+1  . ,n+l 


II'  I x p"'  ' V 


These  equations  are  already  In  the  tridlagonal  form, discussed  In 
section  2.3,wh1ch  allows  for  their  rapid  solution.  This  can  be  displayed 
more  directly  by  combining  terms  In  (3.5)  above  to  yield 


•-x  4-l,j"”x  = 0 


(3.6a) 


S ^1 J-1  * ”y  '^i!]  ^ ^1 J+1  ^ 


(3.6b) 
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where 


•■x  “ ”x  “ ^ '^x*  '^x  ' i^x^l.j  “T'^x 

* “ ^i.j  * 1 Vl,j^'l.j+1  " ^1.j-l^  “I'^y^'i.j+l  “ ^'?,j  * ^IJ-l^ 


% “ “ 'i.j  ■*■  i^x“i,j  ^'t+lj  ^t-l.j^  “7‘^x  ^^1+1 J ~^^U3  * ^i-lJ^ 


C»Ai;c=~;D*  5^  ; D * 

X AX  ’ y Ay  ’ X Rax^  ’ y RAy^ 


(3.8) 


The  linearized  values  of  u and  v are  determined  from  the  solution  of 
the  Poisson  equation  for  the  stream  function,  using  the  auxiliary  relations 
(3.3)  which  define  u and  v.  As  In  the  primitive  variable  procedure,  only 
the  steady  solution  Is  desired,  so  the  Poisson  equation  need  not  be  solved 
exactly  at  each  step.  Only  a solution  for  v*,  reasonably  consistent  with 
the  current  solution  for  c.  Is  required  for  the  method  to  proceed.  Hence 
there  Is  no  time  advantage  to  an  Implicit  solution  procedure  for  the  solution 
of  the  Poisson  equation  since  the  exact  solution  (which  couTd  be  obtained  more 
quickly  via  the  use  of  Implicit  methods)  Is  not  required.  For  the  present 
purposes  either  an  SOR  or  explicit  Integration  technique  was  used.  Both  can 
be  written  as 

♦l.j  “ ’’'l.j  4 *1-1  ,j  W ^*1J+1  ^*1J  ♦l.j-r 

- Ax2  (3.9) 


where  a Is  the  relaxation  factor  for  SOR  or  equivalent  to  4A‘i/ax2  If  an 
explicit  scheme  Is  used,  where  at  Is  a false  time. 
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3.3  Solution  Algorithm 

The  solution  procedure  using  the  ADI-SOR  combination  Just  described  is 
as  follows:  (a)  Compute  the  linearized  velocities,  u and  v,  from  the 
last  computed  values  of  ii;  (b)  Integrate  the  vorticity  equation  one  complete 
AOI  step,  subject  to  the  imposed  boundary  conditions  discussed  below  in 
Section  4,  using  the  standard  tridiagonal  solver,  at  each  of  the  two 
intermediate  steps  which  comprise  a complete  ADI  step;  (c)  Find  a 
new  value  of  ip,  consistent  with  the  value  of  c found  in  step  (b),  by 
making  three  to  five  SOR  passes  (or  explicit  time  steps)  through  the  Poisson 
equation  for  ip;  (d)  Continue  this  process  until  a satisfactory  convergence 
is  obtained. 

3.4  Boundary  Conditions 

The  difficulty  with  the  ♦ - c formulation  of  the  Navier  Stokes  equations 
is  that  the  natural  boundary  conditions  we  have  are  on  u and  v,  not  ip  or 
c.  Both  of  these  velocity  conditions  can  be  easily  related  to  conditions  on 
4),  but  neither  one  directly  relates  to  c.  Thus,  the  boundary  conditions  on 
the  velocity  normal  to  the  boundaries  leads  directly  to  a specification  of  ip 
at  all  points  on  the  boundary.  However,  unless  we  know  the  vorticity  exactly, 
as  in  the  case  of  a specified  velocity  profile  at  an  inflow  boundary;  or  can 
specify  its  normal  gradient  equal  to  zero,  as  a result  of  a downstream 
continuation  outflow  boundary,  there  is  no  analytic  specification  of  c at  a 
boundary  which  can  be  made  without  recourse  to  the  equations  themselves  and 
other  boundary  conditions.  But,  to  complete  the  specification  of  the  Navier- 
Stokes  problem  to  be  solved, a condition  for  the  z variable  on  the  boundary 
(either  c or  its  normal  derivative  given)  must  be  specified. 


For  Inflow  boundaries,  as  previously  stated,  given  the  velocity  profile. 

It  can  easily  be  Integrated  tangential  to  the  boundary  to  get  !|»,  and 
differentiated  to  specify  c*  At  outflow  boundaries,  the  standard,  second 
order  accurate,  one  sided,  first  derivative  form  can  be  used  to  specify  both 
di|)/3n  and  ac/8n  equal  to  zero,  where  n denotes  the  normal  to  the  boundary. 
Taking  for  example,  x to  be  the  direction  of  outflow,  the  following  relation 
for  c (or  i|*)  can  be  obtained  which  does  not  destroy  the  tridlagonal  form  of 
the  matrix  to  be  solved. 


'iMAXj  * T <^IMAX-l,j  ”1  'lMAX-2,j 


(3.10) 


This  Is  In  the  general  form  used  In  (2.34)  to  Incorporate  boundary  conditions 
Into  the  tridlagonal  solver,  and  the  Implementation  Is  the  same  as  described 
following  eqn  (2.34)  In  section  2.4. 

For  solid  walls,  we  can  still  specify  by  knowing  that  the  walls  are 
streamlines  and  If  there  Is  any  mass  Injection  (normal  velocity)  through  the 
walls.  However,  since  we  do  not  know  the  normal  derivative  of  the  tangential 
velocity  (although  we  do  know  the  velocity)  the  vorticity  Is  not  known  dlr^tly. 

A condition  for  must  be  made  which  Is  consistent  with  the  known  boundary 
conditions  and  the  governing  equations.  This  can  be  done  In  a variety  of 
ways,  see  Roache  [11]:  here  we  use  the  most  straightforward.  The  unused 
condition  on  the  tangential  velocity  Is  combined  with  the  governing  differential 
equation,  written  In  difference  form  on  the  boundary  In  question,  to  generate 
a relation  for  the  wall  vorticity,  Zw  ss  follows.  The  Poisson  equation  for 
eqn  (3.4),  can  be  written  at  any  boundary.  In  normal  and  tangential  coordinates 
as. 


* w-  ‘ ^ 
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In  most  cases  (all  that  are  considered  here),  Is  a specified  constant, 
so  that  Its  tangential  derivative  along  the  wall  Is  zero.  Writing  the 
remaining  second  derivative  In  difference  form,  centered  at  the  wall  gives 


(3.11) 


In  obvious  notation.  The  value  of  either  or  Is  not  known, 
depending  upon  the  location  of  the  boundary,  I.e.  the  end  or  beginning  of 
the  domain.  However,  we  doi  have  a relation  which  Involves  this  external 
point,  namely  the  stream  function  definition  of  the  tangential  velocity,  q^. 


In  difference  form  this  Is 


dtv# 


- 1* 

“zsfr 


(3.12) 


Solving  eqn  (3.12)  for  either  or  and  placing  It  Into  eqn  (3.11) 
yields  the  boundary  relation  for  either 

♦w“*- 


S,“ln 


An 


(3.13a) 


at  the  end  of  the  domain,  or 


S»  An 


at  the  beginning  of  the  domain. 


Aii^ 


(3.13b) 


These  are  only  first  order  accurate  conditions  for  but  have  been 
shown  to  be  sufficiently  accurate  In  most  cases.  These  relations  for  c 
are  decoupled  from  the  current  solution  for  ip  and  lead  to  the  sequential 
solution  algorithm  given  In  Section  2.3. 
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4.0  BIHARMONIC  FORMULATION 


4.1  Equations  and  the  Newton-Chord  Method 

Eliminating  the  vortlclty  c In  the  vorticity-stream  function  formula 
tion  yields  what  we  call  the  blharmonlc  formulation: 


Ir  = ffM 


where 


(4.1) 


B.  ^ V 3y  ''  ax  3x  ^ ay 

For  the  steady  state  we  have  simply 

■ 0 


(4.2) 


Newton's  method  for  solving  such  a nonlinear  problem  assumes  given  a current 
approximation  to  the  solution,  say  i|i  « F,  and  then  we  seek  a correction,  say 
such  that  = F + ^ satisfies  (4.2)  when  second  and  higher-order  terms 
in  4)  are  dropped.  This  procedure  gives  the  linear  problem  for  4: 

t(F)*  = -X(F)  (4.3) 

where  the  linear  operator  L(F)  Is  defined  by: 

i(F)4  (4.4) 

Of  course  the  Indicated  procedure  can  be  Iterated  (replacing  F by  F -»■  41, 
which  Is  Newton's  method)  and  terminated  when,  say  ||  4 ||<  10*‘^/2  If  d 
digits  are  required.  This  method  has  many  virtues; not  the  least  of  which  Is 
rapid  convergence.  Specifically  the  error  at  the  next  application  Is 
essentially  ||  ||^  If  F Is  sufficiently  near  a solution. 

The  main  difficulty  with  Newton's  method  Is  that  at  each  Iterate  a new 
linear  problem  (4.3)  - (4.4)  must  be  solved.  However,  the  so-called  modified 


E 


Newton  method  or  Chord  method  can  be  used  to  circumvent  this.  Specifically, 
the  linear  operator  Is  "frozen"  at  say  jC(Fq)  the  first  Iterate,  and 

we  solve  In  place  of  (4.3)  the  system 

Zq*  « -X(F)  (4.5) 

Now  only  the  right-hand  side  need  be  updated  at  each  Iteration  and  the  main 
work  In  solving  (4.5)  need  be  done  only  once.  It  Is  easy  to  show  that  this 
procedure,  although  not  quadratically  convergent,  does  converge  quite  fast 
If  F^  Is  reasonably  near  a solution.  Techniques  have  been  developed  to 
Insure  that  this  combined  Newton-Chord  Iteration,  or  slight  modifications  of 
It  do  converge  and  that  no  more  than  two  full  Newton  Iterates  need  be  done  I 
In  the  formulation  and  discussion  above,  we  have  neglected  the  boundary 
conditions  which  must  be  Imposed  along  with  (4.3)  or  (4.5).  However,  we  will 
consider  only  linear  boundary  conditions,  which  will  usually  be  i|;  and  Its 
normal  derivative  prescribed  (but  may  also  Involve  the  second  and  perhaps 
higher  normal  derivatives).  Then  It  Is  always  possible  to  arrange  that  each 
Iterate  satisfies  the  boundary  conditions  exactly,  so  that  the  correction, 
merely  has  to  satisfy  the  corresponding  homogeneous  boundary  conditions. 

4.2  Difference  Approximations 

Vfe  cover  the  relevant  rectangular  domain  In  the  x,y  plane  by  a grid 
system  defined  by 

Xi-x,  + (1-l)h  (1-1 m),  yj-y,  + (J-l)k  (j  • 1 n) 

(4.6) 

where  h,k  are  the  grid  sizes  In  the  x,y  directions  (see  Fig.  2).  vie  use 

standard  second-order  centered  difference  approximations.  Thus  It  Is  con- 

2 2 — — 

venlent  to  Introduce  the  central  difference  operators  6|^,  6^,  6^  and 

1 
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Figure  2.  A typical  computational  molecule  at  a point  where  an  exterior 
unstored  grid  value  Is  used.  The  dotted  lines  Indicate  the 
grid  values  connected  by  the  two  boundary  equations  (4.27). 
Stored  and  unstored  grid  values  In  the  calculation  are  Indicated 
by  e and  o. 
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also,  to  simplify  the  derivation,  the  corresponding  shift  operators  S,  T, 
In  terms  of  which: 


3 S - 21  + S"\ 


Tk  5 S - S"\ 


s T - 21  + T"^ 


?,3T-r 


(4.7) 


Then  to  approximate  (4.3)  at  each  Interior  grid  point,  I.e.  for  1-2 

m-1,  j - 2,  ...,  n-1,  we  make  the  following  operator  replacements: 

lx  * k^h  = k ly  • 2F*k  = 2ir  ^ 

7^  : (^  ^h  * T = k ^ * k 

h^  k h h K X 

. L (s2  -4S  . 61  -«■’  + S-^')  + V 

(4.8) 

+ _2  (ST  - 2S  - 2T  + ST"'  + 41  + S"H  - 2S"^  - 2T"^  + S H ^ 
h^k^ 

W = <7*6  * ^ ‘h  ' ♦ 2S ’ - s 2) 

+ J-y  (ST  - 2S  + ST"'  - TS"’  + 2S"'  - S'’t"’  ) 
2hk^ 

+ -l-y  (ST  - 2T  + S'H  - ST"^  + 2T"^  - S"H‘^ ) 

2kh^ 

Corresponding  to  the  operators  Jf  and  Ht  these  replacements  generate  a 
nonlinear  difference  operator  K and  a linear  difference  operator  ^ 
given  by 
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“jfc  Vij  ^7  *h^h  * 7 

^ TO  Vij^^  *h*k  ^ J?  *^k^^j 


(4.9) 


(4.10) 


^ TO  *h*k  ^ ^ j*h  4hk  ^ 4*h^ 

"TO  ^h^h  ^ J?  ^ TO  ^h^j  ^ ®k^h^ 


In  terms  of  these  the  difference  equations  for  the  (|>^j  read 


L(F,jU,j  . -N(F,j)  (4.n) 

Using  the  explicit  expressions  In  (4.8)  for  the  operators  occurring  In  (4.10), 
we  can  readily  find  the  coefficients  multiplying  the  4i..  ..  In  (4.11), 

1+Mij+V 

where  n,v  « 0,  +1,  +2  with  |y|  + |vl  <^2.  Denoting  these  coefficients  by 
(see  Fig.  2),  so  that  (4.11)  reads 


2 “^i.ii.j.v  * 

i|  + |v|^2 


(4.12) 


and  defining  for  convenience 


• _4  ^l.^\  2 __4  /1.1\  1 

*2  — T i~7*  -Tit  ■ — jrT»  «a  — 7 (-y  + -y) t a-  * — j-  , 

^ Rh^  h^  k^  ^ ^ Rlf^  Ir^  5 B 


R k^  h*^  k' 


”1  ■ :;x  '’2 " :;i3r  <'h*kf>  * rjr  <'k*k''>’  ‘2  ■ jk  * ■7>'k^' 


h‘  k‘ 


1 T 


(4.13) 
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we  ftnd  that: 


I 

I 


I 


■ ag  + a4  + 2(a^  + a^). 

a°2  * a,  + b, , of,  - -a2  ” *>2  - Cg.  «.]  * *3  + " ^3  • 
®2  * ®1  " h * “1  * ‘*2  ^2*  • *3  “ **3  " *^3  • 

n -1  1 

Cl  * a^  “*  b»* « o **  ^a^  ^ b*  ^ c*  • o » * a<i  ^ ^ > 

“0  5 5 0 4 4 4 -1  3 3 3 * 

*0  " *5  * bg,  o^  ■ -a^  — b^  — c^,  a|  • *3  ” bj  + Cj 


(4.14) 


With  this  notation  we  also  find  that 

“ ®l^h''ij  * *3®h®k^j  * ” Vh^j  * ^4®k^j 


(4.15) 


To  complete  the  system  of  equations  for  the  we  must  adjoin  to  the 
system  (4.12)  the  equations  corresponding  to  the  boundary  conditions.  These 
will  be  considered  in  the  next  section. 


4.3  Boundary  Conditions 

Normally  two  boundary  conditions  will  be  imposed  along  each  of  the  four 
sides  of  the  basic  rectangular  domain  defined  by  the  four  grid  lines  with 
i = 1,  m,  j ■ 1,  n.  We  will  discuss  several  relevant  pairs  of  conditions, 
but  for  simplicity  will  restrict  attention  to  one  boundary  only,  the  one 
along  the  grid  line  i ■ m.  The  corresponding  details  for  the  other  three 
sides  will  follow  in  an  obvious  way.  Certain  differences  have  to  be  observed 
between  the  i ■ 1,  m and  the  j ■ 1,  n boundaries  when  we  consider  the 
solution  of  the  difference  equations,  but  these  need  not  concern  us  at 
present. 

The  most  coimon  pair  of  boundary  conditions  is  for  ^ and  its  normal 
derivative  d^i/dx  to  be  prescribed  functions  of  the  tangential  variable  y. 
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If  the  Initial  guess  Is  set  up  with  the  correct  boundary  values  of  the 
first  condition  Is  simply  represented  by  the  equation 


I 

I 

I 

i 


L 


*m,j  " ® (J  ■ ....  n)  (4.16) 

If  the  same  standard  difference  approximation  as  In  Section  4.2  is  used,  the 
second  condition  can  be  represented  by 


where  q * 2h(34»/3x)  or  In  Newton  form  with  ^ » F + ♦: 

■♦m-l,j  Wu  “ “'"iiH-l.j  ^m-l,J  ^ 

Now  this  condition  Involves  a grid  value  outside  of  the  rectangle  of  stored 
values,  which  must  be  eliminated  by  using  the  appropriate  stream-function 
equation  written  for  one  point  in  from  the  boundary.  Here  this  reads 

(4.19) 

where  only  the  relevant  terms  have  been  specifically  written  down.  The 
remaining  terms  on  the  right-hand  side  do  not  Involve  j*  Elimination 

Vi,j 

...  + (aj  + 02^Vl,j  “iVj  " ^m-l,J^ 


“o*m-l,j  “l*m,J  “2*m+l,J  * ■®2''nH-l,j 


Thus  the  second  boundary  condition  can  be  accounted  for  by  modifying  the 
Newton  form  of  the  stream-function  equation  for  grid  points  one  line  In  from 
the  boundary.  This  modification  amounts  to  replacing  a®  by  a®  ■ a®  -t- 
and  evaluating  the  right-hand  side  In  the  standard  way  by  (4.15),  except  that 
wherever  would  occur,  we  use  In  Its  place  P + P,„_l j* 

If  the  second  normal  derivative  Is  given  Instead  of  the  first,  (4.17) 
will  be  replaced  by 
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1.21] 


where  now  q » h^{aS/3x^)|^a^»  and  the  corresponding  Newton  form  reads 


^*m,j  * Wl.j  * '^m+l,j  ^^ra,j  * ** 


(4.22: 


Again  we  must  eliminate  41^1  j now  between  (4.19)  and  (4.22);  the 
modified  equation  corresponding  to  (4.20)  reads 

...  + (o®  — a«)^  + (a?  + 2a»)^-  j ■ -ao(q  + 2F  4 “ F-  1 4)  + ...  (4.23 

0 2'^nj_l  j 1 Z'^ro.j  2 ^ m,J  m-i,j 

So  again  the  boundary  condition  1s  accounted  for  by  modifying  the  Newton  form 
of  the  stream-function  equation  at  the  appropriate  grid  point  and  by  evaluat- 
ing the  right-hand  side  using  (4.21)  with  ij;  replaced  by  F wherever 
occurs. 

It  may  also  be  relevant  to  Impose  two  boundary  conditions  Involving 


'^m+1  y example,  a satisfactory  way  of  dealing  with  the  downstream 


2 2 

boundary  condition  In  developing  channel  flow  Is  to  Impose  » d t|»/3x  = 


at  stations  suitably  far  downstream.  If  these  stations  are  on  the  grid  line 
1 - m,  we  have 

♦m+1  .j  ” '*’m-l  ,j  “ Vl  ,J  “ ^♦m.j  Vl  ,j  “ ° 

Eliminating  j between  the  two  equations,  we  can  replace  the  second  by 

♦nl.j  * Vl,j 

which  In  Newton  form  reads 

(4.26 


'*m-l,j  *m,j  " '^m-l,j  '^m,j 


■ F. 


- F. 


We  can  now  merely  add  this  equation  In  place  of  (4.16)  and  account  for  the 
first  of  (4.24)  by  (4.18)  with  q « 0. 

It  may  be  noted  that  the  pairs  of  conditions  so  far  considered  can  all 


be  Included  as  special  cases  If  we  assume  that  the  pairs  of  boundary  differ- 
ence equations  can  be  put  In  the  form 


It  will  be  seen  In  the  next  section  that  equations  of  this  form  do  not  dis- 
rupt the  band  structure  of  the  Newton  matrix;  In  fact,  we  will  see  that  the 
first  of  the  pair  could  also  Involve  *^^^3  j without  disrupting 

the  structure,  but  this  Is  not  necessary  for  the  cases  we  consider. 

Note  that  the  first  equation  of  (4.27)  Is  applied  for  j = 2,  ...,  n-1, 
I.e.  at  n-2  points,  and  hence,  totalling  the  corresponding  equations  for 
the  other  sides  and  also  the  Newton  equations  for  the  stream  function,  we 
have  MN-4  equations  without  those  corresponding  to  the  second  equation 
of  (4.27).  The  latter  must,  therefore,  total  2M  + 2N-4,  since  the  total 
number  of  unknowns  Is  MN  + 2(M-2)  + 2(N-2).  Thus,  at  each  corner  we  must 
choose  just  one  of  the  associated  sides  to  contribute  a boundary  condition 
corresponding  to  the  second  of  (4.27).  In  most  applications  It  will  not  matter 
which,  so  we  choose  arbitrarily  to  associate  the  comers  with  the  North  and 
South  sides  (j  = 1,  n)  (see  Fig.  2). 

The  notation  In  (4.27)  Is  transferred  Immediately  to  the  boundary  condi- 
tions for  the  other  four  sides,  a multiplying  the  value  at  the  boundary 
point,  b and  c that  at  the  first  point  Inwards,  and  d that  at  the  next 
point  Inwards.  Thus  with  superscript  W,  S,  E,  N referring  to  the  West, 
South,  East,  North  sides,  we  write  the  set  of  boundary  conditions  as 


•^0,3  - ‘I  ^ * ''i,j  ^ 

^,0  ’ * ^*^,2* 


Fi  j * r'^  c'^Fo  i + d'^F 


2,j 


3,j 


>■  * <=  fl.2  " 


1.1 

^m,j 


t4.28) 


The  values  of  q,  a,  b,  r,  c,  d for  any  particular  case  can  be  set 
out  In  tabular  form  as  In  section  5. 
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4,4  Matrix  Structure  and  Solution  by  Band  Solver 


k 

h 

f 


1 


There  are  various  ways  available  for  the  solution  of  the  linear  system 
(4.11).  Several  block  elimination  schemes  can  be  formulated  and  also  Iterative 
schemes.  Including  those  of  AOI  type.  Since  one  would  like  to  take  advantage 
of  the  Chord  or  special  Newton  method,  a scheme  In  which  an  LU  decomposition 
Is  performed  would  seem  to  be  preferable.  Block  elimination  schemes  can  be 
designed  to  do  this,  but  because  of  Its  potentially  greater  stability,  espec- 
ially for  large  Reynolds  numbers,  we  choose  to  explore  the  practical  applica- 
tion of  a standard  band  solver  with  partial  pivoting  for  stability  (this  may 
be  especially  Important  since  we  hope  to  work  entirely  In  IBM  single  precision). 

The  standard  way  of  organizing  the  system  (4.11)  as  a banded  system  Is 
to  order  the  unknowns  by  rows,  I.e.  In  the  order  4i  ]*•••»  1>  4<|  2*  *** 

The  matrix  elements  for  the  system  (4.11)  plus  the  boundary  equations  of 

the  general  form  (4.28)  has  the  appearance  shown  In  F1g.  3.  The  super-  I 

i 

scripts  N,  S,  E,  W Indicate  to  which  of  the  North,  South,  East,  West  sides  \ 

the  second  boundary  conditions  belong.  For  clarity  the  suffixes  attached  to  I 

the  a's  are  shown  on  one  row  only,  and  the  modifications  to  the  a's  due  j 

to  eliminating  the  first  boundary  equations  are  only  Indicated  by  position. 

Thus,  those  needing  modification  at  the  W and  E boundaries  are  shown  with 
a bar,  those  at  the  N and  S by  a hat,  and  those  at  both  by  a bar  and  hat. 

For  use  In  the  band  solver,  the  diagonals  have  to  be  stored  as  the  columns 
of  a rectangular  array,  B say,  so  they  are  numbered  from  1 to  = 4m  + 1. 

The  main  diagonals  of  the  individual  blocks  have  salient  locations,  which  are 
denoted  by  £.2  and  These  are  shown  In  a sample  row  In  Fig.  3,  which 
also  shows  an  equation  count  k In  relation  to  the  successive  blocks. 

It  Is  Important  that  the  rows  are  scaled  so  that  their  norms  are  of 
fairly  uniform  size  before  the  LU  subroutine  Is  called.  The  natural  scaling 
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of  equations  (4.11)  Is  good  enough  since  the  diagonal  element  Is  constant 
and  Is  likely  to  be  the  dominant  element  In  each  row.  However,  the  rows  with 
the  unit  diagonal  elements  are  normally  much  smaller,  so  scaling  Is  definitely 
necessary,  especially  for  single  precision,  as  tests  have  Indicated,  and  It 
Is  convenient  to  multiply  them  through  by  a^.  In  some  problems  these  unit 
diagonal  elements  do  not  arise  If  the  known  boundary  values  of  are  not 
Included  as  unknowns. 

The  major  storage  requirements  are  for  the  matrix  B,  whose  dimensions 
are  (mn,  and  a matrix  C,  say,  for  storing  the  L of  the  LU  decompo- 
sition, whose  dimensions  are  (mn,  £2)*  where  I2  * ZiiH-l.  In  addition,  we 
need  three  vectors  of  dimension  mn  for  F,  the  right-hand  sides  and  the  Inter- 
change permutation.  This  yields  mn  (Sm+S)  as  an  estimate  for  the  storage 
requirements.  Note  that  the  solve  procedure  SOL  after  the  decomposition 
needs  both  B and  C,  so  that  the  Chord  method  needs  the  same  amount  of 
store  as  the  full  Newton  Iteration. 

Operation  counts  yield  the  following  asymptotic  estimates  for  the  cost  of 
each  stage  of  an  Iteration  step: 

LU  » BM^N,  SOL  » 6M^N,  COEFFS  * 15MN 

Here  LU  Is  the  procedure  for  factoring  L(F)  Into  the  LU  form,  and  COEFFS 
denotes  the  procedure  for  evaluating  the  Newton  matrix  and  right-hand  sides. 
Thus,  neglecting  the  cost  of  COEFFS,  the  ratio  of  a full  Newton  to  a Chord 
step  Is  about  4N/3  and  as  M Increases  It  would  seem  more  efficient  to 
switch  to  the  chord  method.  However,  this  depends  on  the  rate  of  convergence 
of  the  linear  Iteration.  This,  In  turn,  depends  on  the  closeness  of  the 
Initial  guess  from  which  the  Jacobian  matrix  Is  calculated.  If  a sequence  of 
Reynolds  number  cases  are  performed,  with  or  without  extrapolation  from  one 
to  the  next,  this  will  In  turn  depend  on  the  Reynolds  number  spacing.  Other 


strategies  can  be  envisaged*  for  exaaple*  two  full  Neirtons  could  be  used 
each  step  with  chord  Iterates.  If  necessary,  either  between  or  after.  It 
Is  not  really  clear  what  an  optlMua  strategy  slight  be  since  this  Is  rather 
problem-dependent  and  also  depends  on  the  accuracy  required.  For  the 
applications  considered  here  we  choose  a fixed  absolute  accuracy  of  four 
decimal  places  In  the  stream  function  as  a reasonable  objective.  Since 
the  velocities  are  obtained  by  dividing  stream  function  values  by  the  step 
lengths,  this  corresponds  roughly  to  graphical  accuracy  In  the  velocities. 
The  vorti city  can  vary  considerably  In  magnitude,  so  It  would  seem  reasonable 
to  allow  the  accuracy  In  the  vorticity  to  be  defined  Implicitly  by  the 
accuracy  In  the  stream  function. 


5.0  TESTS  AND  COMPARISONS 


A series  of  tests  and  comparisons  of  the  basic  methods  has  been  made 

and  Is  continuing.  Considerable  experience  has  been  gained  In  designing 

Iteration  strategies  for  good  convergence  and  It  has  been  clearly  demon- 

j strated  that  the  techniques  developed  avoid  the  convergence  problems 

associated  with  large  mesh  Reynolds  numbers.  It  should  be  pointed  out, 

however,  that  the  present  restriction  to  uniform  grid  codes  does  put  a 

severe  limitation  on  the  accuracy  that  can  be  obtained  for  large  Reynolds 

numbers.  A subsequent  Investigation  will  examine  the  extent  to  which  this 

situation  can  be  Improved  by  Introducing  a variable  grid  of  an  appropriate 

character.  In  the  present  report  we  concentrate  on  convergence  and  consider 

two  regimes.  The  first  covers  a Reynolds  number  range  for  which  we  may  expect 

2 

accurate  results  by  h -extrapolation.  The  second  deals  with  higher  Reynolds 
numbers  for  which  the  grids  that  we  can  reasonably  use  are  not  fine  enough  for 
2 

h -extrapolation  to  be  of  much  significance,  but  for  which  convergence  Is 
still  possible,  because  we  avoid  the  mesh  Reynolds  number  problem,  and 
results  can  be  obtained  that  give  a reasonable  qualitative  picture.  The 
Introduction  of  a variable  grid  should  give  this  regime  more  qualitative 
significance  and  also  extend  It. 

The  basic  test  problems  are  those  of  the  driven  cavity  and  entrance 
flow  In  a channel.  The  former  problem  has  been  used  as  a test  by  almost 
all  workers  In  this  area  and  thus  Is  essentially  mandatory.  In  addition,  V 

It  Isolates  a basic  feature  which  occurs  In  more  complicated  flows  that  we 

f 

wish  to  consider  later.  The  channel  flow  problem  Is  perhaps  closer  to  some 

i 

of  the  practical  problems  we  wish  to  study,  and  It  also  contains  a free  j 

''downstream”  boundary  on  which  the  development  of  appropriate  "soft" 
boundary  conditions  can  be  studied.  i 
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5.1  Driven  Cavity 

We  consider  the  blharmonlc  formulation  first.  Here  the  boundairy  cond1> 
tions  are  very  simply  dealt  with  since  the  fact  that  u ” v ■ 0 on  all  sides 
except  one,  say  the  North,  where  u » 1 and  v « 0,  Implies  that  the  coef> 
ficlents  In  the  boundary  equations  In  (4.27) are  all  constants  and  mostly  zero. 
In  terms  of  we  have  ® * 0 on  the  East  and  West  sides,  ♦ ■ ipy  ® 0 

on  the  South  side  and  li  - 0,  1*  » 1 on  the  North  side.  In  terms  of  the  net 
function  F^j,  the  conditions  on  the  North  side,  for  example,  read 

SO  that  ■ c*'*  ■ d***  * 0,  q*'*  = 2Ayt  a^  * 0,  b^  * 1.  A similar  treatment 

of  the  normal  derivatives  on  the  other  sides  yields  the  table 


A wide  range  of  Reynolds  numbers  from  25  to  2000  were  tried  on  a very 
coarse  11  x 11  grid  with  a variety  of  Iteration  strategies  to  gain  experi- 
ence wtth.  the  convergence  properties  , Several  strategies  produced  convergence 
for  R > 1000,  but  none  of  those  used  did  so  for  R * 2000.  It  did  not  seem 
worthwhile  exhausting  all  possibilities  In  order  to  see  whether  R * 2000 
was  an  upper  bound  for  this  grid,  but  It  does  seem  likely.  In  any  case, 
although  the  results  for  R > 1000  gave  a qualiutive  picture  that  was  not 
unreasonable,  this  Reynolds  number  Is  well  beyond  that  for  which  one  could 
expect  reasonably  accurate  results  with  a coarse  uniform  grid. 


Several  finer  grids  were  tried  and  since  31  x 31  seemed  to  be  a reason- 
able upper  limit  on  the  IBM  370  In  the  present  context*  most  of  the  tests  were 
performed  on  11  x 11,  21  x 21.  and  31  x 31  grids.  Tests  showed  that  the 
sequence  R > 50,  100,  200  was  a reasonable  range  to  consider  for  the  regime 
In  which  h -extrapolation  could  be  expected  to  give  graphical  accuracy  In  the 
u velocity  on  the  centerline  x * 1/2,  and  the  results  are  graphed  In 
Fig.  4-  The  curves  for  R » 50,  100  are  extrapolated  from  21  x 21  and 
31  X 31  and  seem  adequate  for  graphical  accuracy.  For  R « 200  an  extrapola- 
tion on  all  three  has  been  performed  and  even  this  Is  probably  not  quite  adequate 
for  graphical  accuracy:  the  curve  for  31  x 31  Is  shown  for  comparison.  A 
distinctive  structure  Is  beginning  to  appear  at  this  Reynolds  number.  Results 
for  R > 500  on  the  31  x 31  grid  are  also  shown;  they  Indicate  that  the  sharpen- 
ing up  of  the  structure  Is  being  followed,  but  an  extrapolated  curve  Is  not 
given  since  errors  are  likely  to  be  substantial  to  graphical  accuracy. 


Data  for  various  Iteration  strategies  are  given  In  the  following  table. 


Iterations 

Required 

Convergence 

Factors 

mi 

m 

R • 

Grid  1 

50 

100 

200 

4M/3 

50 

100 

200 

Newton 

Total 

Average 

11 

X 

11 

1 + 

5 

9 

16 

15 

0.26 

0.59 

0.77 

0.010 

0.001 

0.31 

0.010 

21 

X 

21 

1 + 

7 

11 

27 

28 

0.41 

0.64 

0.81 

0.138 

0.003 

0.314 

0.105 

31 

X 

31 

1 + 

7 

12 

33 

41 

0.44 

0.67 

0.85 

0.625 

0.C12 

1.270 

0.423 

31 

X 

31 

1 + 

7 

1 + 2 

1 + 4 

41 

0.44 

0.07 

0.53 

0.625 

0.012 

2.061 

0.687 

31 

X 

31 

E! 

1 + 2 

1 + 5 

41 

0.39 

0.14 

EO 

0.625 

0.012 

2.656 

0.885 

The  number  of  Iterations  Is  given  In  the  form  1^  12»  or  1^  If  l-j  ■ 0, 

where  1^  Is  the  number  of  Newton  Iterations  and  I2  the  number  of  chord  \ 

Iterations.  The  general  strategy  was  to  specify  two  Integers  n-i  and  n2* 

and  to  perform  a maximum  of  n-j  Newtons  for  the  first  Reynolds  number  case 
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and  a maximum  of  02  Newtons  for  the  remaining  ones,  followed  In  a11  cases 
by  chord  Iterations  until  convergence  (max  |4|  < % 10~^).  At  first  It  was 
thought  that  a reasonable  strategy  might  be  to  take  n^  * 2,  02  * 1,  as  In 
the  last  line  of  the  table,  but  evidently  this  would  be  efficient  only  for  a 
much  larger  spacing  In  Reynolds  numbers  than  Is  Justified  by  the  grid  llmlta- 
tatlons.  For  the  sequence  considered  here  It  clearly  Is  most  efficient  to 
evaluate  the  Jacobian  only  once  and  keep  It  fixed  for  three  Reynolds  number 
cases.  It  Is  unlikely  that  one  could  determine  an  optimum  sequence  of 
Reynolds  numbers  In  advance.  It  Is  more  likely  that  one  would  specify  a short 
sequence  of  Reynolds  nuMbers  for  preliminary  Investigation  before  continuing 
to  higher  numbers.  Further,  when  a variable  grid  code  Is  Introduced,  one  may 
very  well  wish  to  change  the  grid  structure  before  continuing  to  higher  values. 
In  such  circumstances.  It  Is  of  Interest  to  consider  the  convergence  from  a 
zero  Initial  guess  to  a higher  Reynolds  number.  The  following  table  gives  some 
data  from  some  early  experiments  with  larger  Reynolds  numbers  (•  denotes  non- 
convergence of  the  chord  Iterations). 

Iterations  Required 

Grid  250  500  1000  2000 


R 


Total  Time 


It  Mould  be  needlessly  expensive  to  make  this  table  more  complete*  but  some 

observations  may  be  made  from  selected  Items.  In  going  straight  to  R ■ 500, 

for  example,  we  see  that  more  Newtons  may  be  needed  as  the  grid  Is  refined; 

and  for  the  sequence  R > 250,  500,  1000  with  the  11  x 11  grid  the  strategy 

n-j  - 2,  n2  * 1 ultimately  falls  for  a Reynolds  number  for  which  a solution 

# 

exists. 

The  comparison  with  the  i>  — i ADI  scheme  of  section  3 Is  strictly  a 
comparison  between  methods  for  solving  the  difference  equations  since  the 
difference  approximations  are  Identical.  What  we  are  comparing  Is  a - c 
Iteration  scheme  In  which  the  c equations  are  solved  by  AOI  and  the  f 
equations  are  solved  by  SOR  with  a blharmonlc  scheme  In  which  the  Iteration 
matrix  Is  factorized  with  a big  band  solver.  The  first  thing  that  can  be 
said  Is  that  the  blharmonlc  scheme  needs  vastly  more  storage  than  the  i*  - c 
scheme,  at  least  mn  (6m  5)  words,  which  for  a 31  x 31  grid  Is  approaching 

2 X 10^,  as  compared  with  little  more  than  6mn,  which  could  be  reduced  still 
further  If  the  convenience  of  storing  u and  v was  relinquished.  Clearly 
one  must  look  for  a compensating  advantage  In  the  blharmonlc  scheme.  The 
boundary  conditions  are  somewhat  more  easily  dealt  with  In  the  blharmonlc 
scheme,  but  this  Is  fairly  marginal.  The  obvious  remaining  characteristic 
Is  speed  or  run  time  to  achieve  a given  accuracy.  This  Is  not  quite  so  easy 
to  compare  as  might  be  expected,  since  the  - t Iterations  can  be  very 
slowly  convergent.  In  which  case  the  same  tolerance  test  max  |a^|  < % 10'^ 
would  not  give  a fair  comparison.  This  Is  because  the  actual  error  can  be 
many  times  larger  than  the  Iterative  change  In  cases  of  very  slow  linear 
convergence.  If,  as  usually  the  case,  the  convergence  Is  ultimately  geometric 
with  a coninon  ratio  less  than  but  close  to  unity,  a comparison  with  the  geo- 
metric series  s^  ■ a + ar  + ...  + ar,  n-1  for  which  |S«  - 11m  S^| 

> |S|^  - S^_i|  X r/(1-r),  suggests  that  we  should  terminate  the  Iterations 
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*(hen  max  |a^|  < % 10"^  x (l-r)/r,  where  r Is  an  estimate  of  the  common 
ratio  of  the  size  of  the  successive  corrections.  Since  the  actual  ratio 
of  successive  corrections  will  not  necessarily  be  monotonic  and  smooth,  we 
use  the  spatial  average  of  the  |ai|i|  and.  In  fact,  store  the  last  eight 
values  of  this  quantity,  say  O]*  >2 *8* 

r *(«5  + Ag  + Ay  + ag)/(ai  + a2  + a^  + a^).  Comparison  of  several  sample 
values  of  the  converged  results  from  each  method  shows  that  they  give  exactly 
the  same  results  for  the  stream  function  to  the  four  decimal  places  required, 
and  also  that.  If  the  original  criterion  of  max  |atp|  < % 10~^  Is  used,  sub- 
stantial differences  In  the  fourth  decimal  do  occur. 

A further  complication  to  a straightforward  comparison  Is  the  choice  of 
pseudo  time  step,  for  experiments  show  that  for  each  grid  and  each  R there 
Is  an  optimum  at  for  which  convergence  Is  reached  In  the  smallest  number  of 
steps.  Unless  we  use  this  optimum  value,  there  Is  no  definite  scheme  which 
we  can  say  Is  typical.  For  the  present  this  has  been  found  by  varying  the 
parameter  y ■ at/ax;  for  example,  y « 4.3  for  R « 100  with  the  31  x 31 
grid  and  thus  yields  a value  of  r « 0.92  for  the  asymptotic  convergence 
factor.  It  seems  that  r Is  less  sensitive  to  y for  y < y^p^  than  for 
y > In  fact,  r Increases  quite  rapidly  for  y > Y^p^*  and  It  may  be 

possible  to  use  this  to  estimate  y^p^  by  Increasing  y as  the  solution 
progresses  until  the  currently  monitored  estimate  of  r exceeds  1.0  and 
then  backing  off  a notch.  For  the  typical  case  mentioned  above  the  time  to 
convergence  (In  84  Iterations)  Is  about  0.125  compared  with  0.625  for 
1 Newton  on  a 31  x 31  grid,  so  there  Is  a fair  amount  of  leeway  that  could 
be  taken  up  In  such  monitorings.  One  could  also  afford  to  store  several 
earlier  Iterations  and  carry  out  some  form  of  extrapolation  such  as  Altken's. 
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Thus  In  the  Reynolds  number  regime  where  h^-extrapo1at1on  Is  accurate, 
the  ADI  scheme  seems  to  be  substantially  faster.  One  presumes  that  the  con> 
vergence  will  get  slower  as  the  Reynolds  number  Increases  and  that  the  problem 
of  finding  an  optimum  y will  worsen.  This  turns  out  to  be  the  case  as  the 
table  below  Indicates,  but  even  for  R ■ 500  the  time  required  Is  still  less 
than  1 Newton  on  the  31  x 31  gi id.  However,  considerable  time  was  used  In 
experimenting  to  find  and  attempts  at  automating  this  have  so  far  not 

been  satisfactory.  Total  time  Is  therefore  becoming  quite  comparable  with  the 
blharmonlc  scheme. 


R 

100 

250 

500 

Optimum  Y 

4.3 

4.5 

4.7 

Convergence  factor 

0.92 

0.96 

0.98 

Iterations  required 

84 

220 

400 

Time  required  (370  mins) 

0.125 

0.322 

0.587 

There  would  seem  to  be  more  point  In  developing  the  variable  grid  codes 
before  trying  to  push  the  Reynolds  number  up  further.  For  variable  grid  codes 
more  time  would  be  spent  In  evaluating  the  residuals,  which  would  tend  to 
Improve  the  speed  ratio  between  the  blharmonlc  and  the  tp  - c schemes,  so 
one  may  expect  the  blharmonlc  scheme  to  be  preferable  with  a variable  grid, 
especially  at  Reynolds  numbers  In  excess  of  500. 


5.2  Channel  Flow  with  Uniform  Parallel  Inlet  Velocity 

Vte  make  use  of  synmetry  and  take  the  channel  wall  to  be  along  x * 0 and 
the  centerline  along  x » %.  Entry  conditions  are  u » 0,  v » 1 along  y * 0 
and  fully  developed  flow  Is  assumed  to  have  been  reached  by  y ■ y,^^. 

Normally  y„„  Is  very  much  bigger  than  >s,  so  we  expect  to  have  to  take  more 
stations  In  the  downstream  direction  than  across  the  channel.  An  appropriate 
y will  depend  on  Reynolds  number  and  In  fact  should  be  proportional  to  R 

nMiX 


for  large  R.  Experiments  showed  that  ■ 3.0  was  reasonable  for  R = 50 

so  for  larger  values  we  took  y„.„  ■ 0.06R.  This  turned  out  to  be  excessive 

max 

for  larger  R.  presumably  because  we  are  not  In  the  asymptotic  regime  at 
R » 50. 

A downstream  boundary  condition  In  which  the  fully  developed  profile  was 
Imposed  generated  oscillatory  behavior  near  the  downstream  boundary,  so  a 
condition  of  parallel  flow  Independent  of  y was  Imposed  Instead.  This 
I avoided  the  oscillations  and  yielded  a downstream  profile  adequately  close  to 

I the  known  parabolic  profile.  Thus,  the  boundary  conditions  read: 

\ X » 0:  u = 0,  V » 0,  X - H:  u ■ 0.  v„  = 0 

y-O:  u-O.  v"1.  “y  ' “ 

or  In  terms  of 


X • 0:  ♦ - 0,  = 0,  X - • 0 

y-0:  ♦■-x.  y-y^,!  ♦y-^yy'O 

Using  standard  central  difference  approximations  and  the  eliminations  mentioned 
In  4.3,  we  obtain  the  following  table 


For  a selection  of  moderate  Reynolds  numbers.  Including  R » 150  for 

[121 

which  we  can  make  comparisons  with  the  results  of  Wang  and  Longwalr  % 
we  obtain  the  following  performance  data,  starting  from  a first  Initial  guess 
given  by  fully-developed  flow. 


1 


1 

1 


5C 


370  Timings  (minutes)  Iterations  Required  Convergence  Factors 


r 


Results  for  the  case  R « 150  on  the  21  x 61  grid  with  * 9.0, 

IiiQX 

Ax  = 0.025,  Ay  - 0.15,  obtained  with  the  biharmonic  program,  were  compared 

[121 

with  the  results  of  Wang  and  Longwalr  \ Satisfactory  agreement  was 
obtained.  A further  series  of  tests  on  a crude  8 x 15  grid  were  made  to 
Investigate  the  convergence  for  larger  Reynolds  number.  The  strategy  defined 
by  n-j  - 2,  n2  = 1 worked  extremely  well  on  the  Reynolds  number  sequence 
R - 25,  50,  100,  200,  400,  800,  1000,  3200,  the  convergence  for  R - 3200 
appearing  to  be  just  as  good  as  for  the  earlier  values. 

A limited  number  of  tests  have  been  made  with  the  41  - c ADI  program 
for  this  case  because  It  soon  became  evident  that  the  b1 harmonic  scheme 
was  clearly  superior  and  a lot  of  pointless  work  would  have  been  carried  out 
In  finding  YQpf  Suffice  It  to  say  that  a typical  particular  case,  R = 50 
on  a 11  X 31  grid,  converged  In  about  230  Iterations  with  y ® 5.0  and  con- 
vergence factor  0.97,  and  took  0.15  mins,  as  compared  with  0.11  mins  for  the 
same  single  case  using  the  blharmonlc  program. 
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